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ABSTRACT 

We model driven, compressible, isothermal, turbulence with Mach numbers ranging 
from the subsonic {M. « 0.65) to the highly supersonic regime {M. « 16). The forcing 
scheme consists both solenoidal (transverse) and compressive (longitudinal) modes 
in equal parts. We find a relation al = blog (1 + b^Al^) between the Mach number 
and the standard deviation of the logarithmic density with b = 0.457 ± 0.007. The 
density spectra follow 'D{k, M.) oc with scaling exponents depending on the 

Mach number. We find ((A4) = aM.^ with a coefficient a that varies slightly with 
resolution, whereas /3 changes systematically. We extrapolate to the limit of infinite 
resolution and find a = —1.91±0.01, /3 = —0.30±0.03. The dependence of the scaling 
exponent on the Mach number implies a fractal dimension D — 2 0.967V1“°'^°. 

We determine how the scaling parameters depend on the wavenumber and find that 
the density spectra are slightly curved. This curvature gets more pronounced with 
increasing Mach number. We propose a physically motivated fitting formula 'D(k) = 
Dgk''^ by using simple scaling arguments. The fit reproduces the spectral behaviour 
down to scales k « 80. The density spectrum follows a single power-law rj = —0.005 ± 
0.01 in the low Mach number regime and the strongest curvature rj = —0.04 ±0.02 for 
the highest Mach number. These values of r] represent a lower limit, as the curvature 
increases with resolution. 

Key words: hydrodynamics, turbulence, method: numerical, method: statistical, 
ISM: structure, ISM: kinematics and dynamics 


1 INTRODUCTION 

Turbulence is an extremely important phenomenon in as- 
trophysical media on nearly all scales. It occurs in ac¬ 


cretion discs (Meschiari 20121, supernova remnants (Inoue 


et al. 20091, star forming molecular clouds (Mac Low & 
Klessen|2004[), the di ffuse phases of the interstellar medium 
I Elmegreen & Scalo 20041, and in the intracluster media 


of clusters of galaxies (Fang et al. 20091. It is known to 


play a dominant role in these environments by influencing 
the morphology, mixing characteristics, statistical properties 
and thermal structure of these gaseous flows. 

The interstellar medium (ISM) consists of complex, 
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chaotic, and filamentary structures, where turbulent mo¬ 


tions interacting with shocks play a major role (Klessen & 


Glover 20141. Radio scattering and scintillation measure¬ 


ments revealed density fluctuations of the ISM consistent 
with a turbulent Kolmogor ov spectrum over a wid e range of 
scales from 10® to 10^®m ( Armstrong et al.]|l995 |. A thor¬ 
ough understanding of these density structures is of great 
importance, as it has a direct influence on quantities like, 
e.g. the star formation rate, star formation efficiency, initial 


mass 

function, and core mass function (Mac Low & Klessen 

2004 

IBallesteros-Paredes et al.||2007[ McKee & Ostriker 

2007 

Hennebelle & Chabrier||2008| 2009 

2011 

1 . Therefore, 


a comprehensive understanding of the mass spectrum ad- 
vected by a turbulent medium is of major interest. 

The theoretical treatment of turbulence by |Kolmogorov| 


(1941a I was derived in the incompressible limit. Hence, the 
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velocity field and its statistical properties were main focus 
in theoretical works, also in the snpersonic, compressible 
context. Only recently theorists started analysing the den- 


Goldreich (19941 and Goidreich & Sridhar (19951 developed 


sity power spectrum in simulations (e.g. 

Padoan et al. 

2004 

Kim & Ryu|2005| [Beresnyak et al.|2005 

[Kowal et al. 

2007 


Kritsuk_et_8jj[200^ Lemaster fc Stone||2009| Schmidt et al.| 

2008[ Fhderrath et al. 112009 1. Kim &: Ryu| ( |2005 1 performed 

a systematic study with a set of simulations finding that the 
density power spectra get shallower with increasing Mach 
number. Federrath et al. (20091 investigated the influence of 


the mode of the forcing routine on the scaling exponent find¬ 
ing that compressive forcing yields a steeper spectrum than 


the solenoidal driven case. Padoan et al. (20041 carried out 


turbulent simulations with weak and strong magnetic field 
finding that the density power spectrum is shallower in the 
strong case. Kowal et al. (20071 confirmed this by studying a 


set of weakly compressible simulations with varying strength 
of the magnetic fields. 

|Saichev &: Woyczynski| ( |1996[ ) developed a model for the 
density field advected by a more general version of the Burg¬ 
ers equation, describing the extreme limit of highly com¬ 
pressible gas or equivalently very large Mach number flows. 


The ’’normal” Burgers equation (Burgers 19391 treats the 


flow as network of interacting shock fronts and therefore ne¬ 
glects the pressure forces (for a detailed description of the 
Burgers equation in the context of the interstellar medium 
see e.g. Klessen fc Glov^2014|. The model of Saichev & 


Woyczynski (19961 permits an analytic treatment of the 


density field of compressible gas with dissipative and disper¬ 
sive effects as well as pressure forces. [Saichev &: Woyczynski| 
(19961 derived analytically that the density spectra follow a 
power law oc fc”"’, with different exponents C, depending on 
the relation between the pressure forces and the magnitude 
of the velocity fluctuations of the medium. They obtained 
a density power spectrum ©(fc) oc k~^ for strong pressure 
forces (i.e. weak shocks) and ©(fc) oc hP in the limit of negli¬ 
gible pressure forces (i.e. strong shocks). Interestingly, these 
are also the results for a density distribution consisting of a 
single discontinuity or an infinitely compressed peak, i.e. a 
step function and a delta function, respectively. Note that 
the power spectrum of a field is determined by its strongest 
singularity. In general a discontinuity in the (g — l)th order 
derivative leads to a power spectrum of the form oc fc“^'* 

( Bee fc Khanin|2007 1. For small values of the Mach number 
Tatsumi & Tokunaga (19741 and Tokuna^ (19761 showed 
that one-dimensional compressible turbulence reduces to the 
solution of two Burgers equations for density fluctuations, 
which implies a density spectrum scaling oc k~^ in the low 
Mach number regime in agreement with the result of jSaich^ 
& Woyczynski (19961. 


We therefore analyse the influence of the rms Mach 
number M on the scaling exponent of the density spectrum 


©(fc, M) oc k‘’^-^'> 


( 1 ) 


guided by the mentioned theories. We note that jGeorge et ah] 
(19841 predicted a pressure power spectrum for the incom¬ 
pressible case scaling like oc fc“^^® using dimensional argu¬ 
ments and confirmed this studying an axis symmetric jet. 


Bayly et al. (19921 analysed the formal asymptotic expan¬ 


sion of the compressible Navier-Stokes equations about a 
uniform state predicting a oc fc“^^® scaling for the density 
power spectrum for weakly compressible flows. [Sridhar 


a theory for mildly compressible MHD turbulence, where the 
density power spectrum developed a Kolmogorov like fc“®/® 
scaling. 

The above mentioned numerical studies focused on the 
impact of different physical processes on the scaling be¬ 
haviour of the density power spectrum. The influence of 
the chosen fitting range, the used regression method, and 
the resolution of the simulation were not or only briefly dis¬ 
cussed. In [Konstandin et ^ ( |2015| ) we demonstrated that 
the resolution and the employed regression method can have 
a major effect on the inferred scaling exponents of the power 
spectrum. We showed that the velocity power spectrum is 
not converged with resolution in the highly supersonic case 
and that the scaling exponents still change dramatically 
with a resolution of 1024^ such that these measurements 
and the inferred trends should not be over interpreted. Be¬ 
side the resolution, the width, and the position of the fit¬ 
ting range were key factors, as the analysed velocity spectra 
were curved instead of following a power law. In the study 
presented here, we therefore analyse the local scaling expo¬ 
nents of the density spectra and measure the variation of the 
power law parameters with the scale fc. We also discuss the 
influence of the resolution in detail, before we describe the 
scaling exponent as a function of the Mach number ({M). 

The paper is organised as follows: Section presents 
the properties of our simulations and Section l^the used 
methods. In Section|4]we show our results and discuss them 
in Section in the context of previous studies. The last 
Section 1^ summarises our findings. 


2 SIMULATIONS 


We use FLASH4 (Fryxell et al. 2000 

Dubey et al. 2008 \ 

with the HLL5R solver (^Waagan et al. 

20111) to numerically 


integrate the continuity equation and the Euler equation 
with a stochastic forcing term F per unit mass on a uniform 
three-dimensional grid. 


dp 

dt 


+ V(p-v) = 0, 


9v , Vp 

_ + (v.V)v = -^+F, 


( 2 ) 


(3) 


where p is the mass density, v the velocity field, and p the 
pressure. We simulate an isothermal medium throughout 
this study such that p = pCs^, with the sound speed Cg. 
This is a reasonable assumption due to efficient cooling pro¬ 
cesses observed in the dense interstellar medium and molec¬ 


ular clouds (Elmegreen & Scalo 20041. To measure the in¬ 


fluence of the resolution on the results, we run simulations 
with 256^ 512^ and 1024^ grid cells. 

We compute the random forcing field F in Fourier space 
with a natural mix of rotational and compressive (solenoidal 
and compressive) modes (Schmidt et al. 2006[ Federrath 
et al.|^10| Konstandin et al. 20121. The forcing only oc¬ 
curs on the large (integral) scales 1 |k| ^ 2, measuring 

fc in units of 2'k/L with L being the box size. The autocor¬ 
relation time-scale of the forcing algorithm is set equal to 
T = L/{csM), where M is the time averaged root mean 
square Mach number. We fix the energy input rate to yield 
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a constant energy flux through the scales. With this forc¬ 
ing scheme we drive rms Mach numbers of At « 1, 5,10,16 
in the statistically steady state of fully developed turbu¬ 
lence for the 256^, 512®, and 1024® resolutions. Additionally, 
we performed 20 simulations with systematically increas¬ 
ing Mach numbers ranging from At = 0.65 up to TVt = 16 
equally spaced with AAt « 0.8 for the lower resolution 256®. 
Therefore, we end up with 32 different simulations in total. 
We start with homogeneously distributed gas at rest and let 
it evolve for > IIT and store the relevant quantities every 
O.IT. 

Figurej^shows the time evolution of the rms Mach num¬ 
ber (left panel) and the standard deviation of the density ap 
(right panel) of the simulations for different resolutions. Af¬ 
ter about two turbulent crossing times the velocity and the 
density fields are in the statistically steady state of fully 
developed turbulence. Both quantities measure the global 
one-point characteristic of the flow and are independent of 
the resolution. We analyse in the following the spectral be¬ 
haviour of the simulations for t 3T. 

Figure [^presents the mass density in a cut trough the 
a;y-plane at 2 = 0 of the simulation 7V(01-1024® (top left), 
A405-1024® (top right), TVd 10-1024® (bottom left), and A416- 
1024® (bottom right) to illustrate the flow pattern in the 
statistically steady state of fully developed turbulence. The 
figures illustrate that the density contrast is increasing with 
the Mach number in general. A closer look reveals that the 
A401 simulation lacks density fluctuations on small scales, 
whereas the higher Mach numbers show density fluctuations 
on all scales. This already indicates that the density spec¬ 
trum is steeper at low Mach numbers, which we will analyse 
in the following in more detail. 


3 METHODS 


The Fourier spectrum of the density field is defined as 

T>(fe, ti, A4)dfc = 47rfe®p(fc, ti, At) • p*{k, ti, At) dfc, (4) 

where p is the Fourier-transformed density field and p* its 
complex conjugate. We use the model described in |Kon-| 
standin et al. (20151 performing a hierarchical Bayesian lin¬ 


ear regression on the logarithm of the power spectra. We 
assume that every time snapshot ti of the density spectrum 
'D{k,ti, At) follows a power law, 

V{k, ti, At) = A{ti, At)fe‘^'**”^VA(fc, At, U ), (5) 

with the amplitudes A{ti, At), the scaling exponents 
C,{ti, At), and the scatter term a^ik, M,ti). This method 
has the advantage that errors and uncertainties are treated 
self consistently, averaging the data is not necessary, as all 
snapshots of the spectra are fitted simultaneously, and it 
provides valid estimates for the fitting parameters, their er¬ 
rors, as well as their time variation. We apply the method 
to extended fitting ranges, as well as small fitting windows 
only containing seven data points [A: — 3, fe -|- 3], which we 
move systematically over the spectrum and interpret the re¬ 
sult as local scaling exponent at the scale k. We refer the 


reader to Konstandin et al. (20151 for a detailed description 


of the Bayesian model, various test on the parameter esti¬ 
mates (like the influence of the fitting range), and a compar¬ 
ison with ordinary linear regression methods applied to the 



Figure 3. Density power spectra of the simulation with 256® 
(dotted), 512® (dashed), and 1024® (solid) resolution and different 
Mach numbers averaged over time t ^ 3T in the statistically 
steady state of fully developed turbulence. 


averaged spectrum as well as the individual power spectra. 
We focus the discussion in the following on the mean group 
slope, its uncertainty and variation with time. The mean 
group slope can be interpreted as a time averaged scaling 
exponent. 


4 RESULTS 


The Parseval theorem 


roo 

(^1+ = / T^{k) dk 

Jo 


( 6 ) 


links the integral of the density spectrum with the variance 
and the squared mean of the mass density probability distri¬ 
bution function. Hence, the first two moments of the density 
distribution describes the area below the density power spec¬ 


trum. We therefore discuss in Section (4.11 the probability 


density function with its moments before we focus on the 
density spectrum and its scaling behaviour in the following 
Sections. Figure [^presents the density power spectra for dif¬ 
ferent Mach numbers and different resolutions. The mildly 
supersonic simulation A401 stands out because of its lower 
amplitude in comparison to the highly supersonic simula¬ 
tions. Figure indicates that the area gain of the density 
spectrum for higher Mach numbers happens through a shal¬ 
lower scaling exponent instead of an increasing amplitude. 
Comparing the resolutions in Figure [^reveals that the spec¬ 
tra with 7V401 are hardly distinguishable for k < 20, whereas 
in the highly supersonic simulations the resolution has an 
influence for k > 5, depending on the Mach number. 


4.1 The probability distribution function (PDF) 
of the mass density 


Passot & Vazquez-Semadeni (19981 propose a heuristic 


model for the density PDF based on the density contrast 
from the shock jump condition in an isothermal medium. 
They assume that the logarithmic density variation As by 
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time t[T] time t[T] 

Figure 1. Time evolution of the rms Mach number (left panel) and the standard deviation of the density dp (right panel) at different 
resolutions 256^ (dotted), 512^ (dashed), and 1024^ (solid). 



Figure 2. The mass density in a cut trough the xy-plane at 2 = 0 of the simulation A401-1024^ (top left), A405-1024^ (top right), 
A410-1024^ (bottom left), and A416-1024^ (bottom right) in the statistically steady state of fully developed turbulence. 
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the shocks can be interpreted as an additive random process 
changing the logarithmic density s = log(p). Arguing with 
the central limit theorem they propose a normal distribution 
for the logarithmic density s and a log-normal distribution 
for the density. The moments of these distributions are con¬ 
nected via 

fj.p = + a^/2) (7) 

for the mean and 

= ffliexp {a^^) - 1) (8) 


for the standard deviation. We choose the parameters of the 
simulations in this way that the box size L = 1 , the total 
volume V' = 1 and the total mass M = 1 such that fip = 1 
implyies 

fis = -oID (9) 

and 


2 / 2 \ -1 
ap = exp ( 0 - 3 ) - 1 . 


( 10 ) 


Hopkins (20131 developed a model for the density dis¬ 


tribution taking mass conservation and intermittent fluctu¬ 
ations into account with a non log-normal shape in order 
to explain the deviations from the relations between the 


moments seen in 3-dimensional numerical simulations (e.g. 

Kowal et al.||2007| 

Kritsuk et al.||2007| Schmidt et al.||2009| 

Price & Federrath 

2010| Federrath et al.||20101 |Konstandin| 


from the log-normal shape can be expressed by a single pa¬ 
rameter T , with T = 0 for a log-normal distribution. 

Figure shows the PDF of s measured in the simula¬ 
tions with 1024® resolution and at different Mach numbers 
together with normal distributions depending only on one 
parameter Gs- We express the mean of the normal distri¬ 
butions with equation The normal distributions are in 
excellent agreement with the measured PDFs also in the 
wings of the distributions. We show the relations ([^ and 
(10 1 between the moments in Figure]^ to quantify the dis¬ 


crepancy between the log-normal assumption and our simu¬ 
lations. Our data show only negligible deviations from both 
relations indicating that our density PDFs are closer to a 


log-normal shape than the ones analysed by Hopkins (20131. 
Therefore, we hnd a T parameter scattering around zero 
(e.g. T = —0.0003 ± 0.0003 in the simulation A415-1024®) 
for all Mach numbers, which we calculated with the relation 
between the moments and the T parameter given in equation 


( 6 ) of Hopkins (20131. The main difference between the sim¬ 


ulations presented here and these in Federrath et al. (2008 1 


and Konstandin et al. (20121 is the decomposition of the 


forcing scheme. We use a here a forcing held, which contains 
both solenoidal (transverse) as well as compressive (longitu¬ 
dinal) modes in equal parts, whereas the above mentioned 
studies focus on the the extreme cases of purely solenoidal 
and purely compressive forcing. 


4.2 The (Ts-A 4 relation 

In the next step we analyse the relation between the stan¬ 
dard deviation of the density distribution and the Mach 
number of the turbulent how. [Passot fc Vazquez-Semadeni| 
(1998 1 concluded from the shock jump condition and the 


CO 


CO 


a. 



-6 -4 -2 0 2 4 

S 


Figure 4. Probability distribution function of s measured in the 
simulations with 1024® resolution and at different Mach numbers. 
The shaded areas indicate the one sigma time variations of the 
PDFs. The red dotted lines correspond to a normal distribution 
with only one parameter, as we expressed the mean value via 
fis = - 0 - 2 / 2 . 


central limit theorem that the standard deviation of the log¬ 
arithmic density Gs ~ bA4 with the proportionality constant 
b. This leads to Gp ~ Gs ~ bA4 when Taylor expanding 
equation (10 1 to hrst order. Federrath et al. (20081 analysed 
the inhuence of diherent forcing helds on the relation 

0-3 = log(l -I- b^A4^). (11) 


They assume a linear relation between the standard devia¬ 
tion of the density field and the Mach number ([Padoan et al.| 


19971 and use equation (10 1 to express the standard devia¬ 


tion of the logarithmic density. They found 0 ^ b 1 de¬ 
pending on the decomposition of the forcing with b = 1/3 for 
purely solenoidal forcing and b = 1 for purely compressive. 
This parameter b can be interpreted as the ratio between the 
compressive Mach number over the total rms Mach number 
of the flow ( Konstandin et al.|2012 | 


The longitudinal modes of the velocity held contain collid¬ 
ing and dispersing flows, which cause density Huctuations, 
whereas the transverse modes have no inhuence. 

The jump condition p{ti+i) = indicates that 

we have to express the PDF at the time ti+\ conditioned on 
the past time p(p(ti+i)| piti)). A stochastic processes is said 
to be Markovian, if the system only depends on the state 
of the previous time step, but not on those before. As the 
shock jump condition suggests this behaviour, we will use 
the Fokker-Planck transport equation for random variables, 

^p(s; t) = -^[Bi(s; t)p{s\ f))] -b ^^[^ 2 ( 5 ; t)p(s; t)] . 

(13) 

with the drift parameter Bi{s, t) and the dihusion coeffi¬ 
cient 52 (s, t) to determine the steady state of the PDF of s 
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0.5 


1 1.5 


2.5 


Figure 5. Relation between the moments of the density field 
assuming a log-normal distribution. Top panel: Mean of the log¬ 
arithmic density s versus its squared standard deviation together 
with equation Bottom panel: Squared standard deviations of 
the density p and the logarithmic density s together with equation 




| |Risken||l996| [Baudoin||2014[ ). An equivalent description of 
the process can be achieved by expressing the random vari¬ 
able itself in a stochastic differential equation instead of the 
time evolution of its PDF 


ds = Bi{s, t)dt + B 2 {s, t)dW , 


(14) 


with the random process dW. Following the idea of |Pas-| 


sot & Vazquez-Semadeni (19981 and Federrath et al. (20081 
we assume Bi — 0 and B 2 ~ log (1-I-b'^TVf'^) in equation 
(141, which describes the density fluctuations as random 
process caused by the shocks. With this ansatz the Fokker- 
Planck equation has no stationary distribution. Instead we 
get with the initial values to and so and the initial condition 


p{s; to|so; to) = (5(s - So) 


p{s; t|so; to) 


1 

y/2Tr(t — to)(Ts 


exp 


2(j|(t -to)) ’ 


(15) 


which is the dispersion relation of a Brownian motion. This 
is expected as the model does not contain the redistribution 
of the mass according to the hydrodynamical equations. In 
the next step we assume that the pressure difference be¬ 
tween the ambient/average pressure pp and the pressure 
p of the different positions of the medium relaxes the gas 
compressed by the shocks. We express the pressure term 
s ~ cRs — fj-a) and make the ansatz 


Bi{s{t)) = - 1 / tq(s - ps) bI = l/rp log (1 -I- b^A4^) , 

(16) 

with the efficiency parameter I/tq and l/rp for the differ¬ 
ent processes, which will be determined later. Equation (141 
reads then 


ds{t) = -{s - ps)/Tadt+ V(log (1 -I- h'^M'^)/Tp)dW , (17) 
which is a Ornstein-Uhlenbeck process (|Pope|2000[ |Baudoiii] 


20141 describing the non equilibrium evolution of the density 


held. 


The stationary solution with ansatz JB can be calcu¬ 
lated with the Fokker-Planck equation (131, which gives 


p(s, t|so, to) =N{ps, Ps) , 


(18) 


a normal distribution with mean ps and squared standard 


deviation nf = 
the theories of 

T'a/2r;3 log (1-I-b^AI^) in agreement with 

Passot & Vazquez-Semadeni (19981 and 

Federrath et al. 

(|2008[) besides the prefactor Tal2Tp. The 


timescale Ta associated with the pressure term can be inter¬ 
preted as the time for redistributing shocked and diluted gas 
and Tp is the timescale associated with the shocks occurring 
during the time dt. We express the timescale of the pressure 
term with the dynamical timescale = L/2M. Whereas we 
estimate the shock frequency with two times the compressive 
Mach number rp = L/dMc as only the longitudinal part of 
the velocity field contributes to advecting flows, however it 
counts twice for opposed shock fronts. Therefore we end up 
with Tq/ 2 t /3 = b such that 


crj = b log (1 -I- h^M^). 


(19) 


Figurej^presents as a function of the Mach number M for 
the simulations at different resolutions. Also shown are the 
relations proposed by [Passot &; Vazquez-Semadeni| (|1998[) 


and Federrath et al. (20081, both fitting the data for M < d, 
but significantly overestimate erf in the high Mach num¬ 
ber regime. We measure with these models b = 0.24 ± 0.01 
{(js = bA4) and b = 0.27 ± 0.06 (erf = log (1 + b^AI^)) 
limiting the fitting range to A4 <4. These measurements 
are slightly smaller than b = 1/3 proposed for purely 
solenoidal forcing. In contrast the fit of relation (191 is in 


agreement with the data for all Mach numbers. We mea¬ 
sure b = 0,472 ± 0.002, 0.459 ± 0.005, 0,457 ± 0.007 for 
256^, 512^, 1024^ resolution respectively. The new proposed 
model could also explain the measurements of |Konstandin| 
|et al.| ( [?012[ ), who found strong deviations from equation (111 
for purely solenoidal forcing (b « 1/3), whereas no discrep¬ 
ancies where found for purely compressive forcing (b ~ 1). 
We will analyse the influence of the forcing field on this re¬ 
lation in future work to test this hypothesis. 
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Figure 6. Squared standard deviation of the mass density cr^ as 
a function of the Mach number. We show the relation between 
the Mach number and the standard deviation of the logarithmic 
density (red dashed line), = log (1 + b^Af^) (red 

solid line) and the new proposed formula = blog (1 + b^Af^) 
(black lines). The fits of the data with 256®, 512®, 1024® resolu¬ 
tion are solid, dashed, dotted, respectively. 


4.3 The local scaling exponents 

Figure shows the local scaling exponents of the density 
power spectra as a function of the scale for different Mach 
numbers. Recall that every point corresponds to a fit over a 
range of Ak — 6 such that the first points are the results of a 
fit fc £ [4 : 10] centred at kc = 7. The lines connect the mean 
population slope, the thick error bar corresponds to its 1 a 
uncertainty, and the thin error bar states the 1 cr variation 
with time. The dashed lines indicate the results of the lower 
resolution runs (512®) and the solid lines are for the higher 
resolution (1024®). Note that simulations with resolution of 
512® and 1024® the wave numbers above k > 16 and k > 32 


are known to be dominated by numerical effects (Kritsuk 
et al.|[^07| [Federrath et al.||2010[ [Konstandin et al.||2012[ 


Bertram et al. 2015 Konstandin et al. 20151, respectively. 
Therefore, we stop showing the scaling exponents for the 
512® resolution simulations for fc > 28 to keep the plot read¬ 
able and focus the interpretation in the following discussion 
on wave numbers smaller than fc < 32. 

All local scaling exponents in Figure]^ are nearly con¬ 
stant and only increase slightly with the scale fc. This cur¬ 
vature increases with higher Mach numbers, such that the 
difference between the large scale k — 7 and the intermediate 
scale fc ~ 25 local scaling exponent is ~ 0.3 for the A416- 
1024 simulation. The low Mach number run A401 has a small 
range (fc G [7 : 18]) in which both resolutions have compa¬ 
rable and nearly constant local scaling exponents. Whereas, 
the higher Mach number runs show a resolution dependence 
of the local scaling exponents already on scales close to the 
forcing routine, where the scaling exponents of the higher 
resolution simulations are systematically smaller (« 0.2) 
than their lower resolution counterparts. 


Mltd - - 10243 



Figure 7. The local scaling exponents of the density power spec¬ 
tra as a function of the wave number k at 512® (dashed) and 1024® 
(solid) resolution and for different Mach numbers. The local scal¬ 
ing exponents are the results of a Bayesian regression applied to 
systematically moved fitting windows (fc — 3, fc -1 3]. The thick er¬ 
ror bar indicates the uncertainty of the group scaling exponent 
and the thin error bar states the variation of the scaling exponent 
with time (see [Konstandin et al.|2015| for more details). 


Our measurements are in agreement with the results of 
Saichev & Woyczynski (19961 predicting a spectrum oc k~^ 


for weak shocks and oc fc^ for infinitly strong shocks. For a 
given resolution, all the curves of the local scaling exponent 
for flows with higher M lie systematically above the ones for 
lower Mach numbers, with values of about —2 for A4 = 1. In 
the limit of weak shocks the density profile is composed of 
sawtooth or step functions. Contrary in the limit of strong 
shocks the density profile is peaky with few spots contain¬ 
ing most of the mass (i.e. delta functions, see Figure [^. 
Following the ideas of Kim & Ryu {20051 our results can 
be interpreted as follows. The density jump in an isother¬ 
mal shock is proportional to the square of the Mach number 
yielding very large density contrasts. Therefore, conserving 
the total mass, the gas gets more concentrated in shocks 
with increasing Mach number causing a more peaky den¬ 
sity distribution. This leads to a shallower density spectrum 
for the higher Mach numbers. Increasing the resolution has 
the opposite effect. If a simulation with high Mach num¬ 
ber is poorly resolved most of the mass is contained in few 
grid cells creating a peaky density distribution with shallow 
scaling exponent. With increasing resolution substructures 
get refined and mass gets distributed over more grid cells 
causing a steeper density spectrum. 


4.4 The influence of the Mach number and the 
resolution 

In the next step we quantify the influence of the resolution. 
Therefore, we fit the power spectrum over an extended range 
starting at the largest scale, which is not directly influenced 
by the forcing routine (fc = 4), and containing a number of 
points that doubles with the resolution. These ranges are 
fc G [4 : 10], fe G [4 : 17] and fc G [4 : 31] containing 7, 
14, 28 data points for the 256®, 512®, and 1024® simulation. 
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Figure 8. Scaling exponents as a function of the Mach number, 
measured in k G [4 : 10] for the simulations with 256^ (cyan), 
G [4 : 17] for 512^ (light-green), and /c G [4 : 31] for 1024^ 
(dark-green) resolution. The thick error bar indicates the uncer¬ 
tainty of the group scaling exponent and the thin error bar states 
the variation of the scaling exponent with time. The results of 
different authors are shown with grey symbols, as stated in the 
legend. Additionally, the power law fits to the 256^ simulation 
(dotted), 512^ (dotted), and 1024^ (solid) are also shown with 
the fitting values stated in the legend. 


respectively, which are comparable to those widely used (e.g. 
Kritsuk et al.|2d07l I Federrath et al.|20ldj |Konstandin et ap 


2012 [ Bertram et al.|20T5 1. We follow the standard approach 

in the literature and fit a power law to the spectrum. As 
we see in Figure]^ the spectrum is actually slightly curved, 
and we note that this introduces additional uncertainty that 
we need to include in our Bayesian analysis. Therefore, we 
analyse the influence of increasing the uncertainty estimate 
artificially on our results as well as varying the fitting range 
in the Appendix [A] 

Figure shows the scaling exponents as a function of 
the Mach number for the simulations with 256^ (cyan), 512® 
(light-green), and 1024® (dark-green) resolution. Addition¬ 
ally, we added the results of Kim & Ryu ( 2005|) (grey c ircles), 
Kritsuk et al. ( |2007| ) (grey square) and Kowal et al. (20071 
(grey diamonds). The data-point of [Kritsuk et aL[ (20071 
(1024®, compressive driven) is in ag r eemen t with our results, 
whereas the data of Kim&^^R^ (20051 (512®, solenoidal 


driven) and Kowal et al.| (20071 (256'^, solenoidal driven. 


weak magnetic field) are systematically shallower for the 
highly supersonic cases. This is caused by the different forc¬ 
ing routines and confirms the finding that solenoidal forcing 
yields shallower density power spectra than mixed or com¬ 
pressive driven ones (Federrath et al.| 20091. Another rea¬ 
son is the weak magnetic field in the simulations of |KowaT| 
et al. (20071, which is known to flatten the spectra further 
] Padoan et al.|2004 1. 

To describe the functionality of the scaling exponent 
with the Mach number, we perform a Bayesian power law 
fit with two parameters C(^) = aM^ on the results with 
different resolutions. We assume that a, /3, and the error 
on the measured scaliug exponents are normally distributed. 
The result of the regressiou is shown as solid (1024®), dashed 


(512®), and dotted (256®) black lines and parameters are 
listed in the legend of Figure With our value /? < 0 the 
model is in ag reement with the theory of|Saichev fc Woy-| 


czynski (1996 1 as it converges towards the scaling V oc k'^ for 


the strong shock regime A4 —>■ oo. ft also recovers the weak 
shock regime oc k~^ at the_Maclinum^_r7^=^_(;^£/(a)^_^ 


and t he scaling D oc k ^/® (George et al.[l984 Bayly et al. 


19921 at A4 = (-7/3a)®/'®. 


To get an estimate of this functionality in the limit of 
infinite resolution, we assume a model with three parame¬ 
ters, 


C(A 1 , n) = , ( 20 ) 

where n is a factor corresponding to the resolution (2” 256)® 
of the simulation. This model uses the assumption that the 
influence of the resolution on the measurement of the scal¬ 
ing exponent halves by doubling the resolution, which is in 
agreement with our individual fits. We perform a Bayesian 
model fitting the three parameters of the scaling exponents 
of all resolutions simultaneously with the result 


a = -1.91±0.01, 0 = -0.70±0.01, lj = 0.20T0.01. (21) 


This is in agreement with our individual fits shown in Figure 
The sum converges for n —>■ oo towards 2 such that we 
get in the limit of infinite resolution P = ft + 2u) 

((M) = (-1.91 ± 0.01) A4'° ®“° °® . ( 22 ) 


This is a remarkable result, as we can confirm for the first 
time that the trend of shallower slopes of the density power 
spectrum with increasing Mach number is independent of 
the resolution. 


5 DISCUSSION 

5.1 The fractal dimension 


In analogy to the hierarchical structure of the velocity, 
characteristic for incompressible turbulence theory, |von| 
|Weizsacker| ( |195T| introduced a hierarchy of clouds. He pro¬ 
posed a theory describing the density distribution of molec¬ 
ular clouds 


pL/ — l — 1 / 


(23) 


where p,, is the density of a cloud at the level of the hier¬ 
archy v, £„ is the size of the cloud at this level, 7 reflects 
the degree of compression, and / is the volume filling factor. 
He assumes a self similar behaviour of the density field such 
that every cloud contains a certain number of smaller clouds 
and so on, yielding density distributions described by equa¬ 
tion (|23l|. In this picture 7 is zero or one for no or isotropic 


compression, respectively. Fleck (1996 1 extended the work 


of von Weizsacker (19511 and proposed a relation between 


the scaling of the density and the fractal dimension D, 


p{£) oc e(f)f“®®' oc e{£)£ 


(24) 


where we added the unit step function 0{£), as p{£) is only 
defined for positive scales £. The Fourier transformation of 
equation (241 gives 

p(fe) = r(l - 37)fc"^^"®®'\ 


(25) 
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for the magnitude, neglecting the phase, and with F being 
the Gamma function. Using the definition of the density 
spectrum equation 0 and equation 0 we get 


-(1-37) = ^C(X), 

(26) 

which leads to 


1 1 ..;3 

7=3 + gaA 4 , 

(27) 

and finally results in 


D = 2- ]^aMY 

(28) 


Note that equation (28 1 is in agreement with the fractal 
dimension defined in Stutzki et al. (19981, who derived a 
relation between the scaling exponent of the density spec¬ 
trum and the fractal drift exponent (Hurst exponent) based 
on the theory of fractal images. They get D = E + 1 — H, 
with the dimensionality E = 2 of the fractal surface and 
the Hurst exponent H = (<( — 2) /2. In the case of fractal 
images, D can be interpreted as fractal box coverage dimen¬ 
sion. With equation ( |28[ ) we get a fractal dimension of D = 2 
and 7 = 1/3 in the limit M —>■ oo . N ote that this is a spe¬ 
cial case, as for 7 —>■ 1/3 equation (24 1 is 0(£)£~^ = S{£) the 
Dirac Delta function, which Fourier transforms to constant 
magnitude oc fc® with a zero phase instead of equation (251, 
as r(0) is not defined. 

At low Mach numbers the average flow velocity is 
smaller than the sound speed such that it can not produce 
significant overdensities anymore. To describe this transi¬ 
tion, we define a critical Mach number Adcrit below which 
we assume a constant density in equation ( |24[ |. The Fourier 
transformation in equations (25 \ and 0 give a scaling ex¬ 
ponent for the density spectrum of C,{Mcrit) = —2 in agree¬ 
ment with the theory of |Saichev fc WoyczynskT| ( |1996[) . As¬ 
suming 7 = 0 and D = 3 in equations and (|28|, we 
obtain 


Merit ^ ( ~ ) ^ 0 '®® ■ 


(29) 


using the measurements of equation (22 1 . For Mach numbers 


M < Merit it follows that < —2 and 7 becomes negative. 
Therefore, 7 can not be interpreted as compression parame¬ 


ter anymore and from equation (241 we see that the density 


fluctuations are confined to the small scales. This regime 
is dominated by sound waves which have a steep spectrum 
with (( = —7/3 ( [George et al.|19M f The proposed range for 
the fractal dimension is also in agreement with observations 


(e.g. Elmegreen fc Falgarone|1996 Sanchez et al.|2()07 and 


references therein) and simulations (Federrath et al. 2009 
Konstandin et al.|2012 1 suggesting an overall fractal dimen¬ 
sion of interstellar clouds in the range D « 2.0 — 2.7. 


5.2 The curvature of the density power spectrum 

The second conclusion concerns the density spectrum itself. 
The discussion in section has shown that the scaling ex¬ 
ponents of the density spectrum are not only changing in 
between the simulations with different Mach numbers (Fig¬ 
ure 1^ but also in each simulation with the scale (Figure 
[^. Gonsidering a driving force producing a turbulent flow 



Figure 9. Same as Figure]^ additionally showing the fit of equa¬ 
tion ( |30| l as dotted lines. The fitting parameters of the scaling ex¬ 
ponents are collected in the legend. Including the new curvature 
parameter 7 to the model improved the fit significantly although 
7 is in most cases very small. 


with the Mach number A4o = Vo/cs at the scale £q and as¬ 
suming a power law scaling of the velocity fluctuations, 
there is a scale £s at which the velocity fluctuations are of 
the order of the sound speed Cs- This scale is called the 
sonic scale (Vazquez-Semadeni et al. 20031. According to 
our analysis in Section |4.4| we expect a scaling exponent of 
the density spectrum of ^ = aM^ above the sonic scale. 
Increasing the Mach number of the box to A4' > Mo yields 
V' /ca = {£o/£'sY > {£o/£aY = Vo/Ca after reaching a statis¬ 
tically steady state, resulting in a smaller sonic scale £a > £'a 
and a shallower density spectrum < Y ■ However, we can 
find a scale £' given by {£'l£'sY = {£o/£aY ~ Mo- This im¬ 
plies that the turbulent flow with higher Mach number M' 
contains a sub volume of size £! with the Mach number A4o, 
the same dynamical range relative to the sonic scale, and 
therefore the same scaling exponent / as the original box, 
while Y = a{M')^ > ( for the new box of size £o- This 
example illustrates that we have to include a dependence on 
the scale in our model of the scaling exponent, which we will 
do in the following paragraph. 

Guided by this picture and the variation of the local 
scaling exponents with the scale, we make a power law 
ansatz for the scaling exponent C,{M, k) = YM)k'^■ We ne¬ 
glect this variation with the scale in Section 13 which can 
be justified as we investigate in equations ( | 20 [ ) and \22\ 
the scaling behaviour of the simulations only in a limited 
range close to the forcing regime. Using t he s caling expo¬ 
nent C{M) = aM^ established in Section 


4.4 


and add the 

scaling factor kP to account for a curved spectrum results in 


V{k) = VoiM)k° 


(30) 


The fit in fe € [4 : 31] of the density spectra with 1024® 
is shown in Figure 0 and the resulting fitting parameter 
for the scaling exponents are collected in its legend. We per¬ 
formed a hierarchical Bayesian model fitting all spectra with 
different Mach numbers at once such that a and /3 do not 
depend on the Mach number and are consistent with the 
measurements in Section 14.41 The additional factor kP im- 
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Figure 10. The energy injection rate e as a function of the Mach 
number A4 with a power law fit. The resulting fitting parameters 
are given in the figure. 
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6 SUMMARY 

We analyse the properties of turbulence using a suite of 
three-dimensional numerical simulations which are contin¬ 
uously driven on the largest scales. The forcing scheme con¬ 
sists both solenoidal (transverse) and compressive (longitu¬ 
dinal) modes in equal parts. We model driven, compress¬ 
ible, isothermal, turbulence with rms Mach numbers rang¬ 
ing from the subsonic to the highly supersonic regime. We 
find the relation — blog (1 -|- between the Mach 

number and the standard deviation of the density distribu¬ 
tion, which improves the £t significantly. We derive this re¬ 
lation with the shock jump condition and the Fokker-Planck 
equation (Section |4.1[ ). We find b = 0.457 ± 0.007 with the 
new proposed formula describing the mixture of compres¬ 
sive and solenoidal modes of the velocity field, which is in 
agreement with our driving scheme. By employing a hier¬ 
archical Bayesian fitting method, we estimate the param¬ 
eters describing the scaling relation of the density power 
spectrum. The density power spectra follow power laws, 
V oc with a scaling exponent depending on the 

Mach number (Section|4.4|) in agreement with the theory of 
Saichev Sz Woyczynskif j 19961. We find that C(.^) — ctM^ 
with a scattering slightly with resolution, whereas /3 gets 
systematically shallower. We model that effect and extrap¬ 
olate to the limit of infinite resolution (equation |20[ | to find 
C(M) = (-1.91 ± 0.01) >i-0'30±0 03 validate this result 
by testing the influence of varying position and width of the 
fitting range, as well as the uncertainty of measured scaling 
exponents of the density spectrum on the inferred parame¬ 
ters (Appendix [A|). 

The dependence of the scaling exponent on the aver¬ 
age Mach number of the density spectrum implies a depen¬ 
dence of the fractal dimension on the Mach number (Sec¬ 
tion 5.11. In the proposed model the fractal dimension is 
D — 2— 1/2C(A4) = 2-|-0.96Af“°'30. The fractal dimension 
is D = 2 in the strong shock regime and D = 3 in the incom¬ 
pressible limit, which is reached at the critical Mach number 
Merit ~ 0.86. This is in agreement with the observations of 


proved the fit significantly such that the fit and the density 
spectra are hardly distinguishable for k < 80. The htting 
parameter r} describes the deviation from a pure power law, 
which is nearly zero for the low Mach number simulation, 
such that in this case the density power spectrum is well 
described by a single power law. For the high Mach num¬ 
ber cases the parameter rj increases up to —0.04 ± 0.02. We 
interpret the measurements of r] as lower limits of the cur¬ 
vature, as the density spectra get still steeper at small wave 
numbers k k, 7 with resolution, whereas their scaling ex¬ 
ponents at intermediate scales changes only slightly with 
resolution (Figure]^ or Figure A1 panel in the middle and 
on the right). 

Before we close this section, we want to mention that 
the parametrisation of our simulations with the Mach num¬ 
ber is arbitrary and we chose it as it is used in the derivation 
of [Saichev fc Woyczynskl] (|1996[) and most studies of super¬ 


sonic turbulence (e.g. Kim &: Ryu|2005 Schmidt et al.|2006 


Konstandin et al.|2012 1 . We can change the parametrisation 

to any quantity that has a bijective mapping to the Mach 
number. Figure shows as examples the energy injection 
rate of our forcing routine, which we fitted with e oc M^'^^■ 
Other possibilities are the standard deviation of the den¬ 
sity, which is CTs = blog (1 + ]PM^) (shown in Figureor 
the compressive Mach number as discussed in [Konstandinj 
jeFallpoT^ . Expressing the scaling exponent of the density 
spectrum with the energy injection rate (via its influence 


Elmegreen & Falgarone (1996 and references therein) sug¬ 
gesting an overall fractal dimension of interstellar clouds in 
the range D ~ 2.0 — 2.7. 

We also determine how the parameters depend on the 
wavenumber and quantify the deviation from a pure power 
law by moving the fitting range systematically over the den¬ 


sity spectrum (Section 4.31. This analysis reveals that the 


on the Mach number, see Figure 101 has two advantages. 
First, it measures the total energy flowing through the cas¬ 
cade as the turbulent boxes are in statistical steady state 
and shows explicitly that we do not need the Mach num¬ 
ber as additional parameter to describe a supersonic tur¬ 
bulent flow. Hence, our results for the density spectrum of 
supersonic turbulence suggest that Komogorovs second hy¬ 
pothesis holds, which states that all small scale statistical 
properties are uniquely and universally determined by the 


density power spectra are slightly curved. This curvature 
gets more pronounced with increasing Mach number. The 
density spectra are steeper close to the forcing scale, shal¬ 
low at intermediate scales and again steeper on small scales. 
The height of this bump in the local scaling exponents in¬ 
creases with the Mach number. 

We develop a physically motivated htting formula re¬ 
producing the deviations from a pure power law based 
on the Mach number dependence of the scaling exponent 


scale i and the mean energy dissipation rate e (Kolmogorov 
|1941a|b| . The second advantage is that it offers additional 
interpretations of our results. 


of the density power spectrum (Section 5.2 I. We propose 
ID(fc) = T>ok‘''^^ with C) = aM^ and a new parameter tj de¬ 
scribing the deviation of the spectrum from a pure power law 
with Hxed scaling exponent. This functionality describes all 
density spectra down to wave numbers of fc ~ 80. We mea¬ 
sure 77 = —0.005 ±0.01 in the low Mach number regime 
such that in this case the density power spectrum follows 

































Density structures in supersonic turbulent flows 11 


a single power-law, and we find the strongest curvature 
rj = —0.04 ± 0.02 in the simulation with the highest Mach 
number. 


APPENDIX A: STABILITY OF THE FITTING 
RESULTS 


As stated in the main part, the approach of choosing an 
extended fitting range, which doubles with the resolution 
fits a power law function, where the density spectrum is 
slightly curved. This could lead to wrong uncertainty esti¬ 
mates of S'® the assumption of a power law is only 

approximately true especially for the higher Mach numbers. 
Therefore, we test in the following how stable the presented 
result is against increasing the uncertainty of C(^)' The re¬ 
sults are collected in Table (A11. The first row contains the 
result already presented in equation (211. Using the stan¬ 
dard deviation of the time variation at as the uncertainty 
results in parameter estimates shown in the second, esti¬ 
mating it with 20% of the measured scaling exponent is in 
the third, and with a constant 0.15 is in the fourth column. 
The last two columns contain results, where we analyse the 
influence of the fitting range. Therefore, we use a constant 
width of the fitting range Afe = 6 for all resolutions and 
the error estimate The fifth column contains the result, 
in which the scaling exponents are measured with a fitting 
range [4 : 10] independent of the resolution, whereas the last 
column contains the result fitting the fiat, shallow parts in 
Figure so [4 : 10], [11 : 17], [25 : 31] centred at 7, 14, 28 
for the 256^, 512®, and 1024® resolution, respectively. The 
parameter estimates in Table show that the presented 
results are independent of small variations. They also con¬ 
firm the finding that the scaling exponent of the density 
power spectrum follows Q{M) oc in the limit of in- 

finit resolution. Figure [AT] shows the scaling exponents and 
the fitting curves of the parameters in the columns 2, 5, 6 
of Table [AT] to demonstrate the goodness of the fits. 
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Table Al. Results of the three parameter fit (a, /3, u)) for different estimates of the uncertainty of the local scaling exponents (first four 
columns) and for different fitting ranges (fifth and sixth column). 



fc e [4 

; 10]256, [4 : 
o't 

17]512, [4 : 31]i024 
20% C 

0.15 

fc e [4 : 10],ii 

A: e [4 : 10]256, [H : 17]512, [25 : 31]io24 

a 

-1.91 ±0.01 - 

-1.85 ±0.01 

-1.93 ±0.14 

-1.90 ±0.06 

-1.87 ± 0.01 

-1.88 ±0.01 

d 

-0.70 ±0.01 - 

-0.65 ±0.03 

-0.71 ± 0.06 

-0.69 ±0.06 

-0.75 ±0.01 

-0.60 ±0.01 

LO 

0.20 ±0.01 

0.19 ±0.02 

0.21 ± 0.05 

0.20 ± 0.05 

0.26 ±0.01 

O.ll ± O.Ol 

d 

-0.30 ±0.03 - 

-0.27 ± 0.07 

-0.29 ±0.16 

-0.29 ±0.16 

-0.23 ±0.03 

-0.38 ±0.03 
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Figure Al. Scaling exponents as a function of the Mach number M for the 256® (cyan), 512® (light-green), and 1024® (dark-green) 
resolution simulations. The scaling exponents are measured in varying fitting ranges as indicated in the figure. Additionally, the result 
of three-parameter fit of equation l|20[| collected in Table [aT] (columns 2, 5, 6) are shown as grey lines from left to right. 
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